%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 6 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: due to the use of a finite number of simulations, 
%%% the final results are random and will not correspond 
%%% exactly to those in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Define variables
nobs = 10^2;%We use nobs=10^4 to create Figure 6 in the paper
N = 25;

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%% \nu=4, T=60
nu = 4;
T = 60;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end
fig = figure;
set(groot,'defaultAxesTickLabelInterpreter','latex');  
subplot(3,3,1)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=4$, $T=60$','interpreter','latex');
ylim([0.95 5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 2 3 4 5]);

%%% \nu=4, T=120
nu = 4;
T = 120;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,2)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=4$, $T=120$','interpreter','latex');
ylim([0.95 5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 2 3 4 5]);

%%% \nu=4, T=240
nu = 4;
T = 240;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,3)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=4$, $T=240$','interpreter','latex');
ylim([0.95 5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 2 3 4 5]);
 
%%% \nu=6, T=60
nu = 6;
T = 60;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,4)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=6$, $T=60$','interpreter','latex');
ylim([0.95 2.5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.5 2 2.5]);

%%% \nu=6, T=120
nu = 6;
T = 120;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,5)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=6$, $T=120$','interpreter','latex');
ylim([0.95 2.5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.5 2 2.5]);
 
%%% nu=6, T=240
nu = 6;
T = 240;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,6)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=6$, $T=240$','interpreter','latex');
ylim([0.95 2.5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.5 2 2.5]);

%%% \nu=8, T=60
nu = 8;
T = 60;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,7)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=8$, $T=60$','interpreter','latex');
ylim([0.95 2]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.25 1.5 1.75 2]);

%%% \nu=8, T=120
nu = 8;
T = 120;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,8)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=8$, $T=120$','interpreter','latex');
ylim([0.95 2]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.25 1.5 1.75 2]);


%%% \nu=8, T=240
nu = 8;
T = 240;
rho = N/T;
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
TrueEta = fzero(fun,0.5*(1+nu/(nu-2)));
TrueVarphi = 2*TrueEta^2*(1-rho)/(nu-TrueEta*(nu-2));
EtaHat = zeros(nobs,1);
VarphiHat = zeros(nobs,1);
EtaTilde = zeros(nobs,1);
VarphiTilde = zeros(nobs,1);
for n = 1:nobs
     tau = (nu-2)./chi2rnd(nu,1,T);
     Y = randn(N,T);
     X = (sqrt(tau).*Y)';
     %Student t
     [~,~,nuhat] = mlstudent(X);
     nuhat = min(nuhat,1000);
     y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
     fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
     EtaHat(n,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
     VarphiHat(n,1) = 2*EtaHat(n,1)^2*(1-rho)/(nuhat-EtaHat(n,1)*(nuhat-2));
     %Elliptical
     muhat = mean(X);
     sigmahat = cov(X,1);
     invsigmahat = inv(sigmahat);
     tauhat = sum((X-muhat).^2,2);
     tauhat = tauhat./mean(tauhat);
     fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
     EtaTilde(n,1) = fzero(fun,1);
     VarphiTilde(n,1) = (1-rho)./(1/EtaTilde(n,1)^2-sum((N.*tauhat.^2)./...
                        (T-N+N.*EtaTilde(n,1).*tauhat).^2));
end  
subplot(3,3,9)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
hold on
plot(1,mean(EtaHat),'r*');
plot(2,mean(EtaTilde),'r*');
plot(3,mean(VarphiHat),'r*');
plot(4,mean(VarphiTilde),'r*');
line([0.6,2.4],[TrueEta,TrueEta],'Color','black','LineStyle','--')
line([2.6,4.4],[TrueVarphi,TrueVarphi],'Color','black','LineStyle','--')
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=8$, $T=240$','interpreter','latex');
ylim([0.95 2]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.25 1.5 1.75 2]);
fontsize(fig,16,"points");